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ABSTRACT 

In an attempt to model the accretion onto a neutron star in low-mass X-ray binaries, we present 
two-dimensional hydrodynamical models of the gas flow in close vicinity of the stellar sur- 
face. First we consider a gas pressure dominated case, assuming that the star is non-rotating. 
For the stellar mass we take Mstar = 1.4 x IO^^Mq and for the gas temperature T = 5 x 10^ K. 
Our results are qualitatively different in the case of a realistic neutron star mass and a realistic 
gas temperature of T ^ 10^ K, when the radiation pressure dominates. We show that to get the 
stationary solution in a latter case, the star most probably has to rotate with the considerable 
velocity. 
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1 INTRODUCTION 

Low-mass X-ray binaries (LMXB) are luminous X-ray sources 
composed of a late-type optical companion (mass less than about 
1 solar mass) and a neutron star. About 100 low-mass X-ray bina- 
ries are known now. Neutron stars in such objects are most prob- 
ably old and have a rather weak magnetic field so that an accre- 
tion disk can extend down to the neutron star surface. The rapidly 
rotating gas is decelerating due to viscous friction. The gas then 
spreads over the stellar surface and forms a boundary layer. Here 
most of the energy is emitted in the form of X-rays, whilst its 
amo unt is comparable with the energy generated in the acc retion 
disk jSunvaev & Shakurij|l986l : lsibgatunin & Sunvae\]|l998t) . 

LMXBs can be divided into two different classes. Very lu- 
minous Z-sources (L ~ 0. 1 - IL^m) have relatively soft, two- 
component spectra, while both components can be approximated 
by black bodies wit h color temperatures of about 1 keV and 
2.5 keV, respectively jOilfanov et al.|[2003h . The other less lumi- 
nous sources (L ~ 0.01 - O.OSLedd) are observed in two states: the 
high/soft and low/hard states. The radiation spectra in the soft state 
are similar to those of the Z-sources, while in the hard state they 
are close to the spe ctra of the Galactic black holes in the hard states 
tearret et alj2000h . 

The soft component can be associated with the radiation from 
the accret ion disk, while the hard one is produced in the bound- 
ary layer jMitsuda et Zll 19841 : iGTlfanov et alj|2003h . On the other 
hand, the spectra from the spreading layer depend on the neutron 
star compactness (mass/radius ratio), which determines the gravita- 
tional field at the surface. Therefore, one can get independent con- 
straints on the equation of state of the matter at extreme densities, 
calcu lating the spectra and comparing them with the observational 
data jSuleimanov & Poutan en 2006). 
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A study of the motion of the matter very close to the neutron 
star is also important for understanding the production of quasi- 
periodic oscillations (QPOs) observed in the kHz r ange from a 
number of accreting neutron stars in LMXBs l lvan der Klis 2000|). 
These QPOs may provide direct ways of measuring effects that are 
unique to the strong gravitational-field regime. However, the ques- 
tion about the nature of QPOs is still open, partly because of the 
complexity of hydrodynamical flows in close vicinity of a neutron 
star surface. Thus, detailed studies of the structure of the boundary 
layer plays the key role in understanding the physics in the vicinity 
of a compact object. 

The first model of the boundary layer was proposed bv |Pringl3 

( Il977l) . who considers it as part of the accretion disk. In his model 
the gas is moving parallel to the disk mid-plane and is decelerating 
due to differential rotation and viscous forces. The effective temper- 
ature of the boundary layer appears to be higher than the maximum 
accretion disk temperature, because the size of the BL is smaller 
than that of the disk, whe reas their luminosities are comparable. 
I Popham & NaravatJ ( Il992h identified non-physical aspects of the 
standard a-viscosity prescr iption and developed a mor e physically 
realistic model of viscosity. Nar avan & Pop ham ('1993) proposed a 
self-consistent model of a boundary layer on the surface of a white 
dwarf and accounted for th e hard X-rays observed in cataclysmic 
variables. ' MedvedevI J2004l) studied the radiative accretion onto a 
rapidly spinning neutron star. They considered a quasi-spherical hot 
settling accretion flow and presented an analytical self-similar so- 
lution describing the boundary layer. 

Ilnogamov & SunvaevI ( 1 19991) considered the boundary layer 
as a spreading layer on the surface of the neutron star. They pro- 
posed that matter spirals along the neutron star surface toward the 
poles due to turbulent friction between matter and stellar surface. 
They used a ID approach, averaging all values in the radial direc- 
tion and assuming azimuthal symmetry. To describe the turbulent 
viscosity they used the Prandtl-Karman universal logarithmic de- 
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pendence of the mean velocity on the distance from the stellar sur- 
face and introduced a turbulent viscosity to characterize turbulent 
velocity and turbulent pressure fluctuations. They also assumed that 
the dissipation of rotational kinetic energy causes a strong energy 
release near the bottom of the boundary layer. With these simplifi- 
cations they constructed a semi-analytical model and showed that 
the kinetic energy of the gas is mostly liberated in two belts above 
and below the equator of the neutron star. 

In order to solve t he ID spreading layer problem, 
llnogamov & SunvaevI i 19991) assumed that the initial rotational ve- 
locity in the equatorial plane is very close to Keplerian. However, 
because of the presence of the boundary layer associated with the 
accretion disk, this velocity can significantly deviate from the Kep- 
lerian value. There is no doubt that the behavior of the gas at higher 
latitudes in the spreading layer strongly depends on the conditions 
in the equatorial plane. Therefore it is important to describe the gas 
flow near the equatorial point more accurately. As a first step we 
present here two-dimensional numerical hydrodynamic solutions 
in the neighborhood of the equatorial point. 




Figure 1. Sketch of the boundary layer on the surface of the neutron star. 



2 EQUATIONS AND COORDINATES 

The full non-stationary system of the hydrodynamical equations is 
as follows. The continuity equation is solved in the form 
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where p and U are density and velocity of the gas, and D/Df = 
djdt -I- t7 • V is the advective derivative. The conservation of mo- 
mentum can be written in the form 
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where Fgr = -GM^taif/r^ is the gravitational force (where Mstar 
is a stellar mass, r is a radius-vector), Fys = p"'V • (2vtS) is the 
viscous force, v, is the turbulent viscosity, 'Frad is a radiative flux, 
K is the opacity, and c is the speed of light. The energy equation is 
formulated in terms of specific entropy s. 



Ds 
Df 



: 2y,S' 



1 

P 



-V • Tr. 



(3) 



^6ijV ■ U is the 



where T is the temperature and S = \(Uij + Up) 
trace-less rate of strain tensor and commas denote partial differen 
tiation. In the following we assume that the gas is optically thick 
and can therefore treat radiation in the diffusion approximation, so 
the radiative flux is given by 
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where K = 16(TsbT^ /(3kp) is the radiative conductivity. 

A sketch of the boundary layer on the surface of the neutron 
star is presented in Fig.[T] The gas is accreting in the disk mostly in 
the radial direction R and turns in the Z-direction near the equatorial 
point E. Since the purpose of this paper is to study the gas flow in 
the vicinity of E, we use cylindrical coordinates, and neglect the 
curvature of the stellar surface. We consider a 2D domain limited 
in the radial direction by the surface of the star and the disk zone, 
where the rotational velocity Uf, is close to the Keplerian value Uk 
(see Fig.[2li. 
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Figure 2. Sketch of the calculated domain and of the buffer zones. 



Table 1. Parameters 



quantity 


case 1 


case 2 


^star 


10* cm 


10* cm 


Mstar 


1.4 X W-^ Mo 


1.4 Mo 




3 X 10* K 


10** K 


Tdisk 


1.5 X 10* K 


10^ K 


Pdisk(O) 


0.1 gcm"3 


4 g cm"^ 


V, 


10', ID** cm- s 


10'" cm^ S-' 


a 





1 



i?star is a Stellar radius, Mstar is stellar mass, T^f.^y and Tdisk are temperatures 
of the star surface and the gas in the disk zone, PdiskC) is the gas density at 
the mid-plane in the disk zone, Vt is a turbulent viscosity. 
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3 BOUNDARY CONDITIONS AND BUFFER ZONES 

The boundary conditions represent an integral part of the overall 
solution in that the values on the boundaries both determine and 
depend themselves on the final solution. They must allow the ac- 
cretion onto the stellar surface and the emission of energy. They 
must also simulate the compression of gas near the surface and al- 
low this gas to become part of the surface, while the gas settling 
depends itself on the input parameters of incoming gas and condi- 
tions on the stellar surface. The rotational velocity of the gas in the 
disk part is close to Keplerian, but cannot be exactly Keplerian, be- 
cause otherwise accretion would stop and the boundary layer would 
disappear. On the other hand, the deviation from the Keplerian ve- 
locity is determined by the conditions on the stellar surface. To cope 
with these complications we use so-called buffer zones, which are 
narrow regions just outside the domain (of the size of typically 5 
grid-points; see Fig. where additional terms are added to the 
hydrodynamical equations. This type of approach has proven to be 
useful in earlier simulation s of disk outflows and star-disk coupling 
dvon Rekowski et al.l2003h . We use three different buffer zones that 
are characterized by the three profile functions ({R) for the disk 
buffer zone, A{R) for the star buffer zone, and 7/(Z) for the surface 
buffer zone. In the following we describe the properties of these 
three zones separately. 

In the disk buffer zone, where ((R) = 1, the gas is accelerated 
close to the Keplerian speed due to an additional source term in the 
equation for U^, while Ur adjusts itself to the conditions inside the 
domain. Thus, the (j) and R components of equation ([2j are modified 
by additional terms on their right hand sides. 
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where j denotes the meshpoint in the R direction, Nr is total number 
of the grid points in the i?-direction, and dots indicate the presence 
of terms that where already specified in equation (O- We take t = 
56t, where 6t is the length of the time step. 

In the buffer zone near the star, where A(R) = 1, the radial gas 
velocity goes down to zero at the stellar surface. To describe the 
rotation of the star we introduce the parameter < a < 1, which 
equals zero if the star is non-rotating one, and unity if it rotates with 
the corresponding Keplerian velocity. Thus the and R components 
of equation lO are modified further by the terms 
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where A{R) = 1 in the buffer zone, and zero outside. 

The surface buffer zones will be discussed separately in the 
following two sections, because they have to be treated differently 
for gas and radiation pressure dominated regimes. 

The temperatures on the stellar surface and the disk (left and 
right boundaries of the domain) are fixed by Tstar and T^sk, respec- 
tively, while the gas density is extrapolated on both sides. On the 
lower boundary of the domain (mid-plane of the disk) we use an- 
tisymmetric boundary conditions for the Z-component of the ve- 
locity and symmetric boundary conditions for all other quantities, 
while on the upper domain boundary all quantities are extrapolated. 
The turbulent viscosity v, is assumed to be constant everywhere in 
the domain. 



For all simulations presented here we use the Pencil CodeQ 
which is a high-order finite-difference code (sixth order in space 
and third order in time) for solving the compressible hydrodynamic 
equations CBrandenburg & Dobler, 2002) . 



4 GAS PRESSURE DOMINATED CASE 

As a first test we consider a gas pressure dominated case and choose 
the gas temperature in the disk to be T^a^k = 1-5 x 10^ K. This 
means that the radiation pressure is about two orders of magnitude 
smaller than the gas pressure. Also, we take the stellar mass to be 
Msiar = IO^^Mq so as to balance the gravity force near the surface 
by the gas pressure force. In addition, we assume that the star does 
not rotate (a = 0). We consider two cases with Vt = lO' and 10* cm^ 
s-i. 

Since initially the disk is assumed to be in vertical hydrostatic 
equilibrium, the vertical velocity should be close to zero. In ad- 
dition, we let the gas density p approach a certain vertical profile 
Pdisk(Z), where the value at the disk mid-plane is Pdisk(O) = 0. 1 g 
cm"', and assume that p decreases exponentially with Z. However, 
since the temperature profile results from a thermal balance be- 
tween viscous heating and radiative cooling, the local sound speed 
Cs in Eq. ^ is recalculated at each time step. This allows the ver- 
tical density profile to adjust to the conditions inside the domain. 
Thus, we have in the disk buffer zone 



Dlnp 
Df 

DUz _ 
Dt 



_ Inp - Inpdi 
r 

. - — aR), 

T 



where 

Pdisk(Z) =pd,,sk(0)exp 



and - 



R'cl 



(9) 
(10) 

(11) 



where ((R) = 1 in the disk buffer zone, and zero outside. 

In the surface buffer zone, where 7/(Z) = 1, we assume van- 
ishing first derivatives for all three velocity components t/, {i = 
1,...,3) and for the specific entropy. We also correct the density 
profile to account for the resulting artificial pressure force which 
works against the vertical gravity. Due to this term, the gas flows 
out through the surface boundary rather than coming into the do- 
main at the beginning of the calculation. Thus, we add the terms 

Inp' -lnp^-' -hZV(2//2) 
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We consider two runs with v, = 10^ and 10* cm^ s ' and show 



in Fig. |3] two cross-sections respectively for R 

-i 



1.25 m 



cases. In Fig.|4]we show velocity and density. One can see that the 
accreting gas comes to the stellar surface and turns toward the poles 
of the neutron star. 

We find that the size of the boundary layer, where the rota- 
tional velocity of the gas decreases from the Keplerian value 
down to zero, strongly depends on the value of the turbulent viscos- 
ity. The boundary layer becomes 3.5 times thicker if one increases 

' |http : //www . nordlta . dk/sof tware/pencll- code| 
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Figure 3. Gas velocity, density and temperature as a function of fi and Z for the gas pressure dominated case (Mstar = 1.4 IO^^Mq). The dotted and solid lines 
correspond to Vt = lO' and 10* cm^ s"', respectively. Left panel: fixed Z = 50 m, right panel: fixed R - R^tm = 12.5 m for Vt = lO' and R - Rstai = 20 m for 
Vt = 10** cm- s"' cases. 



Vt by a factor of 10. This is consistent with the classical theory of a 
boundary layer, according to which the thickness of the boundary 
layer is inversely proportional to the squar e root of the Rey nolds 
number, y/Re ~ 1/ (see, for example, Ishih-I Pail[T962l) . The 
increase of viscosity also leads to a growth of the gas temperature, 
resulting from a balance between turbulent viscous friction and ra- 
diative cooling. One can see that the temperature achieves its max- 
imal value in the middle of the boundary layer, where the velocity 
gradient and therefore the heating rate are maximum. 

Note, that the solution for our test case looks similar to the 
spreading layer model for the white dwarf case. We find that the 
Z-component of the gas velocity Vz = 10* cm s"' is very close to 
that obtained bv lPiro & Bildstei] ( |2004|) . Unfortunately, we cannot 
compare other quantities because the values in the spreading layer 
model are averaged along the i?-direction. 



5 RADIATION PRESSURE DOMINATED CASE 

Let us now consider a realistic neutron star mass, Msi^,. = 1.4Mo. 
We find that, in order to balance the gravity force by the gas pres- 
sure gradient near the stellar surface, the gas temperature must at- 
tain a value of 3 x 10'^ K, which is unrealistic. Therefore, in the 
case of a non-rotating star, the only force which can work against 
gravity is the radiation pressure gradient. Indeed, if one takes the 
gas temperature to be T = 3 x 10* K and p > 0.3 g cm"', the radi- 
ation pressure force becomes comparable to the gravitational force 
(i.e. GMstai/R^ g crsBT*t.^J(cp), where ctsb is the Stefan-Boltzmann 
constant). 

Using standard disk theory, we estimate the turbulent viscosity 



near the stellar surface v, = a^Hc^ =^ 10 cm- s" , where ffjisij ^ 
0.01 is a viscosity parameter, H ^ O.Oli^star, and Cs - 10* cm s"'. 
Note, that the radiative viscosity = Aa^^T'^m^l (kc^ p) 10* cm- 
s"' (where is a proton mass) is much smaller than the turbulent 
one, and can be neglected. 

We find that the description of buffer zones must be modi- 
fied in the radiation pressure dominated case. For simplicity we use 
a similar density profile in the disk buffer zone as it was in the 
gas pressure dominated case. However, we now take pdisk(O) = 4 
g cm"' and fix the gas temperature Tjisk to avoid large radiation 
pressure gradients and therefore the generation of large velocities, 
which lead to strongly non- stationary behavior in the buff'er zone 
and eventually to numerical instability. Also, we use a softer con- 
dition for the Z-velocity by assuming vanishing first derivatives. 
Thus, we set 
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In the course of the calculation we find that inside the domain 
cold low-density patches surrounded by denser hot gas appear spo- 
radically. Such patches are dispersed due to motion of the gas from 
the core of the patch outward through the radiation pressure gradi- 
ent. However, if we were to fix the Z and R velocities in the star 
buffer zone to zero, such a patch cannot disappear, while the den- 
sity inside this patch is going to decrease together with the tem- 
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Figure 4. Density and velocity fields in a vicinity of a star in a gas pressure 
dominated case (Msia, = 1.4 IO^-Mq). The domain is limited in the radial 
direction by the surface of a neutron star and in the disk midplane. The 
surface and disk buffer zones are excluded. The viscosity is V| = lO' cm^ 
s"' , the density scale is po = 10"^ g cm"' 



perature, and a numerical instability develops. To avoid this we use 
symmetry conditions relative to the inner boundary of the buffer 
zone R = R(, (see Fig.|2j for Inp, s and Uz- In addition we assume 
Uz = at R = Rq. Thus, we have in the star buffer zone 
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Dt T 
Since the main goal of this paper is to consider the gas flow in 
the vicinity of the equatorial point E, we consider a domain located 
inside the disk with a vertical size smaller than the height where the 
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Figure 5. Temperature and velocity fields in a vicinity of a neutron star in 
a radiation pressure dominated case (Mjtar = 1.4 M©). The domain is the 
same as in Fig.|4]for Vt = lO'" cm^ s"'. 



gas becomes cold and optically thin. To imitate a disk photosphere 
in the surface buffer zone we include an additional cooling term to 
create a vertical temperature gradient and to allow gas to escape 
through the surface boundary. 

First, we attempt to consider a non-rotating star. It turns out 
that in this case the initial distribution of the main quantities (tem- 
perature, density, and velocity) have to be close to the final state, be- 
cause otherwise inhomogeneities in the temperature result in large 
radiation pressure gradients which, in turn, generate large local ve- 
locities. Such velocity perturbations may produce a local decrease 
of density, which will lead to a decrease in radiative cooling, and 
hence to an increase in temperature. The resulting radiation pres- 
sure gradient will decrease the density even further, which leads 
therefore to an instability. (Note that this does not happen in the 
case of a smaller gas temperature because then the gas pressure 
dominates over the radiation pressure.) 
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Figure 6. Gas velocity, density and temperature as a function of R for tlie 
fixed Z = 10 m (dotted curve), Z = 50 m (solid curve) and Z = 110 m 
(dashed curve). The star buffer zone is excluded. 
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Figure 7. The sum of the gravitational, centrifugal and radiation pressure 
forces in the R-direction as a function of i? for values of Z as in Fig.|6] 



Finding suitable initial conditions is a difficult task. It turns out 
that it is easier to consider first a rotating star, and then to decrease 
its rotational velocity down to the necessary value. However, the fi- 
nal rotational velocity still has to be considerable so that centrifugal 
and gravitational forces are of the same order of magnitude. In the 
opposite case, i.e. when the star is almost non-rotating, the gravity 
force has to be balanced by the radiation pressure gradient, which, 
in turn, should be negative near the stellar surface. However, it is 
not clear how to realize this, because the main heating mechanism 
is due to viscous friction which is maximum in the middle of the 
boundary layer rather than at the stellar surface. 

At the current stage we assume that the rotational velocity of 
the star corresponds to the Keplerian velocity at the stellar radius 
(i.e. a = 1). Also, we take a uniform initial distribution of tem- 
perature in the R direction. In that case, at the beginning of the 
calculation, the huge gravitational force near the stellar surface is 
balanced by the centrifugal force. Depending on the temperature 
gradient, the ^-component of the radiation pressure force nearly 
vanishes and the resulting R component of the velocity appears to 
be small. 

While the Z component of the gravitational force is much 
smaller than the R component, its influence on the gas motion in 
the Z direction is crucial and leads to strong flow of gas into the 
domain through the surface boundary and to an accumulation of 
gas near the equatorial point E. To avoid numerical problems, we 
assume that initially the Z component of gravitational force is bal- 
anced by the corresponding component of the radiation pressure 
force. Therefore, the initial distribution of temperature takes the 
form 



r(/?,) = m_i)- 



3GM,t^cp(RdZ, 
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where AR is the size of the mesh in the R direction. We use the 
same initial density distribution as in the gas pressure dominated 
case (see Sec.|4]l. 

In Fig. |5]we present temperature and velocity fields (see also 
the sketch in Fig.[T). Here we show temperature rather than density 
(as was done in Fig.|4ll, because in the radiation pressure dominated 
case the flow of the gas is mostly determined by the radiation pres- 
sure force and hence by the temperature distribution. In addition, 
we take a larger domain size because the thickness of the boundary 
layer now appears to be an order of magnitude larger than that in 
Sec.g] 

The results of the calculation are also presented in Fig. [6] 
where density, temperature, and velocity of the gas are shown as 
functions of R for three diff^erent distances from the mid-plane, 
Z = 10, 50 and 110 m. We find that at higher latitudes the gas 
rotates with a velocity that is comparable to the rotational veloc- 
ity of the star, while in the equatorial plane its rotational velocity is 
smaller than the stellar surface speed. This means that in the equato- 
rial plane the centrifugal force is larger than the gravitational force, 
while at the higher latitudes the centrifugal force is smaller than the 
gravitational force. Such a relation between the main forces would 
result in accretion in the equatorial plane and excretion at higher 
latitudes, provided the gas pressure was much larger than the ra- 
diation pressure. However, since now the radiation pressure domi- 
nates, we obtain the opposite result: the gas accretes only at higher 
latitudes, while in the equatorial plane it is excreting. 

In Fig. [7] we present, as a function of R, the sum of the R com- 
ponents of gravitational, centrifugal, and radiation pressure forces, 
F,o, = Fg,. + Ul,/R + K Tmilc. One can see that at a higher latitudes 
the total force is at some radius negative, F,ot < 0, so the generated 
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radial velocity is negative {Ur < 0), which means accretion, while 
in the equatorial plane Fto, > and Ur> Q,so the gas is excreting. 
The dominant role of the radiation pressure in driving the velocity 
field is also clear from analyzing the temperature in Fig.|5] One can 
see that the temperature decreases outward in the equatorial plane 
and increases outward at Z ^ 110 m near the disk buffer zone. Note 
that along the stellar surface the temperature is almost constant, 
so the Z component of the gravitational force dominates here, and 
causes the gas to sink toward the equatorial plane. 



6 CONCLUSIONS 

We have studied the gas flow in close vicinity of a neutron star in 
a low mass X-ray binary and have assumed that the magnetic field 
is negligible. The main purpose of this work was to investigate the 
flow near the equatorial plane between disk and star, so the curva- 
ture of the stellar surface in the latitudinal direction was neglected 
and cylindrical coordinates were used. 

In the unrealistic, gas pressure dominated case the gas tem- 
perature is r ^ 5 X 10^ K (which is about an order of magnitude 
smaller than the observed value). If the star does not rotate, the 
gravitational force at a radius close to the stellar surface should be 
balanced by the gas pressure force. To have gas pressure and gravi- 
tational forces of the same order of magnitude, the stellar mass was 
chosen to be about two orders of magnitude smaller than the real 
mass of a neutron star. In this case the maximum release of energy 
occurs in the middle of the boundary layer, where the gas velocity 
gradient (and hence the viscous heating) reaches a maximum, while 
the radiation pressure force at the stellar surface is directed inward. 

For a realistic neutron star mass, M^i-^ = 1.4Mo, the gas pres- 
sure gradient at the stellar surface becomes negligible compared 
with the gravitational force. The latter is balanced by the radiation 
pressure, so the gas temperature is about T ^ 10** K. Unlike the gas 
pressure dominated case, the radiation pressure force is directed 
outward rather than inward. Thus, the maximum energy release oc- 
curs directly at the stellar surface. However, it is not clear how to 
realize such a scenario, where the gas is heated by viscous friction 
between the differentially rotating gas layers. 

The picture becomes crucially different if one considers a ro- 
tating neutron star: the gravitational force can now be balanced by 
the centrifugal force. Here we have assumed that the neutron star 
rotates with the Keplerian velocity at the stellar radius. Alterna- 
tively, one might find a solution for smaller rotational velocities by 
gradually decreasing it down to the required value, using the result 
of calculating it for a larger velocity as an initial approximation for 
calculation with the smaller value. 

We find that at higher latitudes the centrifugal force is larger 
than the gravitational force, while at the equatorial plane the gas 
rotates with a velocity that is considerably smaller than the corre- 
sponding Keplerian value. It would be reasonable to assume that 
the gas is accreting near the equator and excreting at higher lati- 
tudes. However, we find the opposite: the accretion occurs only at 
higher latitudes, while in the equatorial plane the gas is excreting. 
This is related to the temperature distribution, and therefore, to the 
radiation pressure force, which is now dominant. We find that near 
the equatorial plane the temperature decreases outward, so the gas 
is pushed away from the surface by radiation pressure. At higher 
latitudes, some distance away from the surface, the temperature de- 
creases inward, resulting in accretion. The circulation of the gas is 
closed by a flow along the stellar surface from high to low latitudes. 



because the temperature is almost constant in this direction and the 
gas flow is controlled only by the tangential component of gravity. 

Finally, one should note that, since we have considered only a 
laminar two-dimensional model of the boundary layer at the neu- 
tron star surface, we have assumed that the turbulent viscosity is 
constant everywhere and that it can be treated as an input parameter. 
Future three-dimensional simulations will allow us to model turbu- 
lent processes more accurately. However, even the results of the 
two-dimensional simulations give us some clues for understanding 
the physical processes near the neutron star surface. These results 
can in principle be used for a more detailed description of the ver- 
tical structure of the boundary layer and for calculating spectra of 
neutron star radiation. Furthermore, the presented results may be 
useful for understanding the nature of quasi-periodic oscillations. 
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